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I. INTRODUCTION 


A. DESCRIPTION OF THE PROBLEM 


The study of dynamic responses in a submersible vehicle using a nonlinear 
analysis is important in determining operating envelopes for the vehicle. Previous studies 
have shown that in straight line motion, the coupling of the sway, yaw, and roll equations 
produce oscillating losses of stability in the system. [Ref.5] Introducing a nonzero pitch 
angle to the vehicles motion will allow us to study the changes to the stability domain for 
a variety of environmental conditions. 

In our work, the linearized equations of motion are studied rants an eigenvalue 
analysis to determine the systems stability through a variety of operational parameters. A 
ronlinear analysis is then conducted to assess the stability characteristics of the resulting 


limit cycles and their impact on the operating characteristics of the vehicle. 


B. OBJECTIVES AND OUTLINE 


In this thesis we expand on previous thesis work which examined the problem of 
stability in straight line motions of a submersible vehicle [Ref. 13]. The primary cause 
for this loss of stability is the coupling es the sway/yaw/roll equations of motion for 
a submersible vehicle. We know that the loss of stability creates stable limit cycles in 
Straight line motion. This work analyzes the effects of introducing a nonzero pitch angle 
to the generic equations of motion in order to determine their effect on the creation of 


limit cycles. 


The model used for this work is a variant of the Swimmer Delivery Model used in 


[Ref. 5] for which a generic set of hydrodynamic and geometric properties are available. 


II. EQUATIONS OF MOTION 


A. COORDINATE SYSTEM 


A moving coordinate system was used for our analysis with the origin fixed at the 
vehicles center of buoyancy. The x-axis is fixed to the longitudinal plane of symmetry for 
the vehicle, the y-axis is positive to starboard, and the z-axis is positive downward. All 


symbols used in the development of the equations of motion are summarized in Table 1. 


B. GENERAL FORM OF THE EQUATIONS OF MOTION 


The equations of motion for a submersible vehicle in the horizontal plane are 


written as follows: 


Sway equation: 


mi + Ur — wp + xg(pq+F)- ye (p? +17)+ 26 (ar - p)I= 
VO eee eee ian 
aoe eV OO Y Ore egal eV eee iect 


(W — B)cos@ sing — 





*nose 2 2 
I. c,, h(x\v + xr)” + Cy b(x)(w xq) | U. G) (1) 


Yaw equation: 


IF +(1,, —Ly)pa-1y(p? -9?)-1,, (pr +4) + L.(ar- p)+ 
mlx, (0+ Ur —wp)- yeU — vr + wa) |= 

N p+ Nona opal 5 qr NGUVceAN vale 
N;U’6,+N,v+N,Up+N,Ur+N,,vq+N,,wp + N,,,wr + 
(x,W —x,B)cos@ sing +(y,W — y,B)sin@+U°N ,,,, — 


i- " c,, A(x\v co xr) + Cp. b(x\(w - xq) eo xdx . (2) 


Roll equation: 
Pp a ~1, lor +1,,(pr—4)-1,,(¢’ -r?)-1,(pq+r)+ 


mly, (w—Ugq + vp)— ZG (v +Ur— wp)|= 


Pe ea I ae Glacial at 
KU’ +K,v+ hee CUT a aD es Wr ot: 


prop 


(yW — y,B)cos@ cos@¢ — (z,.W —z,B)cos@ sing 
(3) 


The rotational velocity equation around the x-axis: 


g=p (4) 


U-, denotes the cross-flow velocity: 


Us (x)= y(v+ ale +(w—xq) (5) 


via 
aa 
M,N) 
m 
p,q?) 
Qy) 
U 
v,w) 
(x,y,z) 
X,¥,Z) 


XG,Y¥cG.zG) coordinates of the center of gravity 
Xp,Vp,Zp)  |coordinates of the center of buoyancy 
Pee fore coordinate of vehicle body 
Da aft coordinate of vehicle body 


Table 1: Nomenclature 





C. SIMPLIFICATIONS 


We must simplify the equations of motion in order to reflect the fact that we are 


analyzing motion about the y-axis. The simplifications that we employ are: 


. Acceleration, w, in the z-direction is zero. 
° Acceleration in the longitudinal direction, u , is zero. 
2 Rotational velocity, g, and acceleration, g, in the y- direction are zero 


D. SIMPLIFIED EQUATIONS OF MOTION 


After applying the above simplifications, the equations of motion become: 


Sway equation: 


mly-+Ur—wp+x-r- Ve (p? noes Ze p|= 
Y,p+Y,r+Y,Uv+Y,,vwt+Y¥, U°S, + Y,v+Y,Upt+Y,Ur+Y,,wpt+Y,, wrt 
(W —B)cos@ sing — 


a Ic b, h(x Xv i xr) +C By b(x}w? v+2r) dx (6) 


\- Ue ( x) 





Yaw equation: 


ieee —I, pr+1,,p+ 
mlx. (v+Ur — wp )— VG (G7 —vr)|= 
N,pt+N,r+N,Uv+N,,vw+N, US, +N,v+N Up+N,Ur+N,,,wp+ 


N,,wr+(xoW —x,B)cos@ sing+(y,W —y,B)sin@+U?N 


prop = 


a Ce h(x\v+xr) +C, b(x)w? (v+sr) 
D: ( \w om (x) (7) 


Roll equation: 


I,,pt+1,,pr—-1,,r° -1,,7+mly, (vp)-z,(v+Ur-wp)|= 


K ,p+K,r+K,Uv+K,,vw+K 


prop 


U*+K,v+K,Up+K,Ur+ 
K,,,wp+K,,,wr+(ygW—y,B)cosé cos¢@—(z,W —z,B)cos@ sin @ (3) 


Roll rate: 


¢=p (9) 





II. LINEAR ANALYSIS 


A. LINEARIZATION 


The simplified equations of motion can be written in the matrix form: 


Ax = Bx +(x) 


where the state vector, x , 1s defined as, 


(10) 


Vv 
c 
< 
p 
Q 
and the state matrices are, 
te mxg-Y, —-mz%,-Y, 0 
~ Bh ie l, —Ne ee 0 
mig — Kon 1 a 0 
0 0 0 J 
and 
ary Dele eu +) W+ im 0 
ie te al XU TN Nw LV Ue ae 0 
UW MgO eh +i w ZW he) tee ay 0 
0 0 ] 0 


The g(x) term remains as given [Ref. 13]. 


The nonlinear terms are linearized around a nominal point: 


X 9 =(venieen Olt =0 


A Taylor series expansion is applied to the nonlinear terms about the nominal point, xo , 


and keeping only the linear components, the equations of motion, written in matrix form, 


become: 
A’x =B’x C1 1 
where, 
AT=A 
and 
YU Y.U-—mU v0) (W -B)cos@ 


IN Ole tix 0 NGO aay (x,W —x,B)cos@ 
KU mz,U+K,U KU .(-z,W+z,B)cos@ 
0 0 | 0 


, 


To assess the dynamic stability of the vehicle, an Eigenvalue analysis is performed in the 


next section. 


B. LOSS OF STABILITY 


An Eigenvalue analysis is used to determine the stability of the linearized system. 
Stability is dependent on the location of the Eigenvalues or the roots of the characteristic 


equation: 


det(B’— AA’ =0) (12) 


in the polynomial form: 


AM +BA’ +CA?+DA+E=0 (13) 


The coefficients equation (13) are given in [Ref. 13, pg. 11]. Using Routh’s criterion we 
can examine the stability of the system. The following conditions must be applied to the 


characteristic equation (13) to ensure all roots have negative real parts: 


BED>ADba=FB- =O (14) 


E>0O (15) 


If E is less than zero, one real root of (13) becomes positive and the system will become 
unstable in a divergent manner [Ref. 9]. This is the case of a directionally unstable ship 
which is well known in the literature [Ref. 3]. If the condition in (14) is not met, it means 
that there is at least one complex conjugate root with real aes and will result in an 


oscillatory response, indicating loss of stability. 


To analyze the limiting case of loss of stability equation (14) must be solved in the 


following form: 


BCD-AD* — EB’ =0 (16) 


The result of this equation will be a curve of zg as a function of xg and will be our locus 


for loss of stability. We can express the coefficients of equation (16) in the form: 


A=A,z,+A,Z, +A; (17) 


B=B,z.+B,z, +B; (18) 


A, Ao, and A3 are of the form given in [Ref. 13, pg. 14]. 


B, and B3 are of the form given in [Ref. 13, pg. 14]. With the addition of a pitch angle, w, 


B> takes the form: 


B, =-m(K,U\I,, -N, )—m(I,, -N; \¥,U) 
+mY, (UN, -Umx, )+mK,(UN, -Umx, ) 
—m(N , +1,, (UY, —Um)+m(N ,U \mx, -Y,) 
—m(K, -1,, XN,U)+mUY, (mx, -N, ) 
—mU(N,, +1, km-Y, )+mUK, (mx, - N,) 


+mz,w(m-Y, XI. —N,)—mz,w(mx, =e \mx —N,) 


C=C,z.+C,z, +C, (19) 


C, and C3 are of the form given in [Ref. 13, pg. 15]. With the addition of a pitch angle, 


w, C> takes the form: 


C, =mU (mx, —N, \Y,U )-mvY, (N,U)—mUK, (N,U) 
+mU(N, +1, X¥,U)—mu (w ,U \m-Y, ) 
+Wim=Y, XI. =i )-W(mx, —Y, Kx. i) 
—m(X ,B-x,W \imx, -Y,)+m(UN, —Umx, \Y,U) 
—m(UY, —Um\N ,U )+m(K,U XUN, —Umx, ) 
—mz,w(m-Y, XUN, —Umx, )—mz,wl¥,U XI,. -N;) 
—mz,w(mx, ae XN ,U)+mz, w(x, -N, XUY, —Um) 
D=D,z,+D, (20) 
D, is of the form given in [Ref. 13, pg. 16]. With the addition of a pitch angle, w, D> 


takes the form: 


D, =UK, (x,B-x,W \m-Y, )-UK, (N,UY,U) 
+UK,(N ,U \Y,U)-(K, -1,, \x,B-x¢W)\Y,U) 
+ K,(x,B-xgW XUY, -Um)-(K,U \x,B—x,W Xx, -Y, ) 
+(K,U \uN. —Umx, \Y,U )-(K,U \UY, -Um XN ,U) 
-(K ,U \Y,U XUN, —Umx, )+(K ,U (UK, -Um\N,U) 


+ mz,w(¥,U XUN, St )-mz,w(UK, —~UmXN,U) | 


The equation for the coefficient, EF, remains unchanged [Ref. 13, pg. 16]. Applying the 
stability criterion, equation (16), and utilizing them in the resulting fifth order 
polynomial, F,, [Ref. 13, pg. 18] we are able to solve F using the MATLAB program in 


Appendix A. The curves show zg as a function of xg and we show results for varying 


pitch angles, w, and varying forward velocities, U. On all of the following graphs the 
pitch angle is varied from 10 degrees to —10 degrees in increments of five degrees. The 


top curve is the 10 degree curve and the bottom curve is the -10 degree curve. 
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Figure 1: ZG vs. XG for U = 2 ft/sec | 
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Figure 2: ZG vs. XG for U = 3.5 ft/sec 
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Figure 3: ZG vs. XG for U = 5 ft/sec 
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Figure 4: ZG vs. XG for U = 6.5 ft/sec 
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Figure 5: ZG vs. XG for U = 8 ft/sec 
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From the results of Figures 1 through 5, the following conclusions can be drawn: 


|. 


In all cases, a sufficiently high metacentric height is required in order to ensure 
vehicle stability. Regions of parameters that fall below the critical curves 
correspond to dynamic instability. 

The critical metacentric height that is required for dynamic stability 1s an increasing 
function of both vehicle speed and trim angle. 

Static stability alone does not necessarily ensure dynamic stability of motion during 
the turn. 

The loss of stability experienced here is a dynamic loss of stability. At the critical 
metacentric height, one pair of complex conjugate eigenvalues possesses a zero real 
part. This is an oscillatory loss of stability, which can not be predicted by considering 


the uncoupled sway/yaw equations of motion alone. 
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IV. NONLINEAR ANALYSIS 


A. INTRODUCTION 


In the previous chapter we performed a linear analysis of the equations of motion, 
observing that with changes to specific parameters (xg, ZG, w) it 1s possible to pass 
through a stable region to a region of loss of stability. 

In previous work [Ref 13] 1t was shown that these bifurcations to periodic 
solutions were all supercritical. This means that limit cycles were produced after the loss 
of stability. By introducing a pitch angle, w, in the nonlinear analysis we will analyze 
whether the bifurcations to periodic solutions will remain supercritical and what changes 


occur to the limit cycles themselves. 


B. THIRD ORDER EXPANSIONS 


Our linearized system was written in the form of equation (11), ignoring the 


nonlinear terms. Including the nonlinear terms changes the form of equation (11) to: . 


A’ =B’x+g(x) (21) 


where: 


wri 
- (x 
BO)=| (x) 


g(x) 


Keeping terms up to third order, the matrix g(x) can be written in the vector form: 


g(x)= g(x)+ 9° (x)+e(z) (22) 


where g)(x) contains the second order nonlinear terms and g(x) contains the third 


order nonlinear terms. The cross flow integrals and the second order nonlinear terms 
remain unchanged with the addition of a pitch angle, w. [Ref. 13, pg.22, 23] However 


the third order nonlinear terms take the form: 


Byaen) 
(3) 
82 (x) 
g(x)= 
gn) 
elec) 
where: 


go) a7) -=(W-B)p" cos@ 
| 
2°) = -]%) -2(x,W x58)" cos @ 
| 
2g?) =< (c,W-zsB)° cos@ 


g;) =0 


The Taylor series expansion yielding the second and third order linear terms of the cross 


flow integrals, and the inverse of the system matrix, (A’)” , remain unchanged. [Ref. 13, 


— —~ 
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pg. 24, 25]. However the B’ matrix has changed as shown in the Linear Analysis section 


of this work. These changes to the 4 x 4 F matrix are: 


Fy Fa Fs Fu 
D D D D 
B= pe). Dp 
D D D D 
a 0 l 0 
where: 
Fi, =4,, Y,U )+a,,(N,U )+a,;(K,U) 
Fis Sms —m)U +a,(N, —mxg U +a,;(K, +z, U 
Fy; aot (Y,U)+ Q1o (NU }+ ays (KU) 


= ae (Ww —B)cos@ eae (x W —x,B)cos@ + a,,(-z,W + z,B)cos@ 
F,, =a,,(¥,U )+a,,(N,U )+a,,(K,U) 

F,, =a,,(¥, -mU +a,,(N, —mx, U +a,,(K, +mz, U 

F,, =a, (Y,U)+a,,(N,,U)+ay(K,,U) 

roa (W —B)cos@ ane (x .W —x,B)cos@ tae = Dey ct z,B)cos@ 
F,, =a,,(V,U )+a3,(N,U )+a,,(K,U) 

F,, =a,,(Y, -mU +a,,(N, —mx, U +a,,(K, +mz, U 

F,, =a;,¥,U;)+a,,(N,U)+a,,(K,U) 


F,, =a,,(W —B)cos@ +a,,(x,W —x,B)cos@ +a,,(-z.W +z,B)cosé@ 
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The remaining elements of the nonlinear and constant terms remain unchanged [Ref . 13, 


Doezone 


C. COORDINATE TRANSFORMATIONS 


To continue our analysis it is necessary to bring our transform our coordinate 
system from state space to a normal coordinate system. This transformation 1s performed 


in the manner given [Ref. 13, pg. 29, 30]. 


D. CENTER MANIFOLD EXPANSIONS 


The center manifold expressions are of the form given [Ref. 13, pg. 30-33]. 


However, the coefficients /;,i=1, 2, 3; j=5, 6, 7 are: 


a 
_ ay 2 2 
L>= YG (m2, +m? ) 


D 


G12 ( 2 ) 
+ D I M3 +1 ,,M3,My, + YVgMymM, 


2 
ay; I ,M3,M,, +1,,my, TMYGgm M3, 
D 2 2 vs 
+y,Wm,, —y,Bm,, > 


Qi; 
Li 6 eG (2m,, m3, +2m,, +my ) 


D 


a, 21 .,M3,M3p ie (m,,my, +Mz)M,, ) 
D 


TYG (m,,m,, +My)M,, ) 


= Gee 


D 


les (m,,my, FM3,7M, )+21,,myMy 
(24) 


+My¢ (m,,m5, + MM, )+ 2(y~W = y,Bynymy 


oe 


ayy 2 2 
1, = D YG (m2, +m, 


G12 ( 2 ) 
a D I M3 +1 ,,.M3.Mzq) + YVgMyM), 


a 2 
os IT, My My, +1).My, +MYgM,M3, +(y,W is y,Bym2, | 


as, ( 2 2 | 
los = D Yg Nz, TM, 


A. 2 ) 
+ D (7,, m2, +1,,M3,M, + YgMymM, 


2 
an [MyM +1,,m3, + MyYegmM Ms, 


L a ygWm;, = y,Bm;, 


oni 
L6 = D Yq (2mm) + 2m», +1, ) 


ae 21 .,M3)M3, tee (m3,mz, + MyM», ) 
p 


D \+ yG (mm, +M,M), ) 


Ay, Lees (m,,m,, TMM», )+21 ,,myMy 


D\+ MYG (m, m3, FMM, )+ 2(ygW a y,B)mymy 


A>, 2 2 
lo 4 = D YG (m2, + m2, 


a2 ( 2 ) 
Bs D I M35 +1, Mz)Myy + YVgM2M 


a3 2 2 | 
= D U7, 259 My, +1,My, +MYgMy Mz, + (y,W a y,Byni, 
a3 2 2 
5 = D DG (m2 +m?) 


a 
32 2 
a D (7,,m2, +1\,M3,M, aie Ygmm,, } 
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(25) 


(26) 


(27) 


(28) 


a 
33 - Wm: : 


a 
l.6= ay YE (2m,,ms, +2m,, +m, ) 


Q3> 


D 


t Vo (m,,m,, +My™M,, ) 


D i+ MY Gg (m,,m5, +M,,Ms3, )+2(ygW ~ y,B)namy 


ae le (m,m, +M3)M4, )+ 21 ,,My,Myy | 
a 
l,,=—y (m? +m; ) 
CAULEY, 22 
D 
G27 m2 eT nym, © ) 
2 D xy N32 TE yz. + YGM M2 


a 
S| mM eer, +MY GMM, +(yoW—y,B)n3,| 


E. AVERAGING 


(29) 


(30) 


(31) 


The procedure for averaging the equations up to the third order 1s the same one 


given in [Ref. 13, pg. 37-40]. The addition of a nonzero pitch angle yields new | 


coefficients J,, =1, 2, 3; j=1, 2, 3, 4. They are: 
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Cp, 





3 2 2 3 
(E,m3 Geen te EI IT 5, +E,m;, ) 


-- - (W —B)cos@m3, 


Nome 


Dy 3 2 2 3 ) 
es (Pmt igre, ita ue . ene, tt: 
as 6y 


=| 
+2 (xqW —x,B)cosOm,, 


— 


a 
+ Ee (zqW —z_B)cos Om, (32) 
oe Cy, SE mt +3E, (m2 m,, a 2m,,m,,m,, ) ] , 
L,=-—+|— +—(W —B)cos03m2 m,, 
a +3E, (m,,3, +2M,)My,M, )+ enue 
Cp, 3E,m,,m, +3E, (mj, my, +2m,,m,.M,, + 3 Ey (m,.m3, +2m,m,M,, 
pe oy emer 
D 
I 2 
+e (oW —x,B)cos 03m: m,, 1 
+2 (2 .W —z,B)cos93m;,m,, (33) 
6D 
Cp, SLB aE +3E, (mj,my +2m,,m,.My, )+ 3E, (m3,m,, +2m MM), 
Ne eae on + 3E,m,m», 
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D 


ie - (W —B)cos @3m,,m2, 


oS 





G l 
‘| a” (E,m3, PS ee ns ay is, }+—(egW -xB)cosn 
Y 


+> (cqW —z,B)cosmi, 3 (34) 


pS “>. (pg 1, + 3E\m, 3E 2 + Eym?, | ‘(Ww B)cos @m; 
2D | 6y 9M), FIL; MM, TIL,M,,M,, +L£L,M), re — B)COS OM, 





ae Co, 3 2 2 3), | 3 
-*2 oy (E,m:, +3E, mM, +3E,m mz, + Em), )+-g (a —x2B)cosema bs 


2 2 2 
Cp, SE m,,m,, +3E, (m2 my, +2m,,M,,M,, )+ 3E, (m,,m?, + 2m ,My,M,, 


Qi oy +3E,m3,m,, 


2 Sy 
+ 7 (W —B)cos03m;.m,, 


2 2 2 
Cy. 3E\m,m,, +3E, (my, 2, +2m,,M,.M», )+ 3E, (m,,m2 +2m,,m,,m,,) 





ay | ©Y |+3E,m2.m,, 


l 
| +g eaW —x5B)cos3m4, ma | 


a3 


+= (2,W 2,8 )cos@3miym C8) 


2 2 
Cp, 3Eym,,M,, +3E, (m3, my, + 2M, ,M,.My) )+ 3E, (m2,m,, +2m,,mym,, ) 


oot 8Y | + 3E,m,,m, 
D 


l3 = 


[+£(W —B)cos 63m, mi 


2 2 2 
Cp, 3E,m,,m), +3£, (m2,m,, +2m,,M,.M, )+ 3E, (m2,m,, is 2m,,m,mM,) 


Go (i +3E 4M, M3, 


| 
\F (x ,W —x,B)cos@3m,,mi, 


L 
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a 
+ (zgW — 9B )cos 03m, ms (37) 


C 1 
lo, = a Gy Eom OE ten +3E,m,m3, +E,m;, )+=(W —B)cos 6m; 


, am (E, m}, +3E,m,m,, oe rtnies +E.) 
ae 
D I 3 
+ (xgW —x_B)cos6ms, 


+3. (2 W—z,B)cos Om, (38) 
6D 
__4u| ©o, 





] 
31 6y (Eym;, + Sein: +3E,m,m,, + Em), )+—(W —B)cos6m;, 


a —- (E,m;, +3E,m,,m,, +3E,m,,m}, +E,m}) 
ees 
D 
l 
+2 (xqW —x_B)cosém,, 


to 


a 
+25 (ceW — 2B )cos Om, (39) 


2 2 2 
Cy, 3Egm,,m,, + 3E, (rj, +2m,,M,,My, +3, (r,.m?, + 2m m»m,,) 


gu oY Th 3E mM» 


D 


ls, a 


+3 (W —B)cos63m:,m,, 


2 2 2 
Cp, 3E\m;,m,, +3E, (m?,m,, +2m,,m,,M, )+ 3E, (m,,m3, +2M,MyM,, 
ay, | SY +3E,m3,m5, 


+ cs (x ,W —x,B)cos@3m:,m,, 


ON 


+22 (co —z,B)cos@m;,m,4, (40) 


2] 


OY aoe (m2, m,, Slee Teves )+ 3E ,M,M5, 


(W — B)cos 03m,,m:, 


os 


—— 


a of oM,,M,> oe (m3,m,, ai 2m,,m,,m, ) 
Cp, 


SEM ip oe (nZ.m,, nF 2m,,m,,M,)) 





6 2 
osEh3 2% y +3E, (m2, m,, + 2m,,m M1, )+ Se gli. 


+ = (%qW ~x_B)cos 63mm 








a a (z,W ~ z,B)cos 63m,,mi, (41) 
C 
“ (E,m}, ai 3E,m,, My, +3E,m,m5, sg E,m3,) | 
] a3) oY 
7 
+iw — B)cos @m;, 
fr) 0 J 
- (E, m3, +3E,m;\m,, SSS gos 3 E,m},) 
yy | OY 
i re (x,W —x,B)cos 0m; 
iG AGW —Xp )cos 42 
a . 
+ on oW — z,B)cos 0m; , oe 


The remaining procedure for determining the equation for the radius of the resulting limit 


cycles is identical to the one described in [Ref. 13, pg. 43, 44], which 1s: 


R=a’eR+ KR? (43) 
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F. LIMIT CYCLE ANALYSIS 


At steady state, R=0, equation (43) becomes: 





0 = R(a’e + KR”) (44) 
Equation (66) has two solutions: 
a (46) 
K 


Equation (45) is the trivial solution and gives no useful information. Equation (46) gives 
a constant amplitude limit cycle, R, in the cartesian coordinate system. This limit cycle 


will exist if the quotient inside the radical sign is positive: 


“aé (47) 
K 





This condition is necessary for R to be areal number. Since a’ is always negative, the 
existence of limit cycles depends on the parameter K. The introduction of a nonzero pitch 
angle does not change the dependence the limit cycle has on the parameter K. Stated, this 
dependence is: 

e If K<0, periodic solutions exist and they are stable. 


e If K>0, periodic solutions exist and they are unstable. 


Zo 


G. RESULTS AND DISCUSSION 


Results for the stability parameter, K, are presented in Figures 6 through 15 
Figures 6 through 10 provide results for K at a constant forward speed, U, and varying 
pitch angles. The pitch angles were varied from positive 10 degrees to negative 10 
degrees in 5 degree increments. In Figures 6 through 10, the bottom curve represents 
solutions for positive 10 degrees and the top curve represents negative 10 degrees. It is 
clear that all values of K are negative indicating they are stable solutions. Notice that 
decreasing pitch angles the solutions for K become less negative, indicating that while 
they are stable, these solutions are less stable than those at the higher pitch angles. In 
Figures 11 through 15, solutions for the stability parameter, K, are again represented, but 
with the pitch angle held constant and the forward speed, U, varied from 2 ft/sec to 8 
ft/sec in 1.5 ft/sec increments. The bottom curve in Figures 11 through 15 represents U = 
2ft/sec. As forward speed increases, the curves representing the stability parameter K 
tend to become more negative, but in a more pronounced fashion than those where U was 


held constant. 
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Figure 6: XG vs. K*GAMMA for U=2 ft/sec 


x10~ U=3.5 ft/sec 





0.4 0.6 0.8 4 ie ae 1.6 1.8 2 
XG 


Figure 7: XG vs. K*GAMMA for U=3.5 ft/sec 
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Figure 8: XG vs. K*GAMMA for U=5 ft/sec 


U=6.5 ft/sec 


0.6 0.8 1 1.2 1.4 1:6 
XG 


Figure 9: XG vs. K*GAMMA for U=6.5 ft/sec 
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Figure 10: XG vs. K*GAMMA for U=8 ft/sec 
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Figure 12: XG vs. K*GAMMA for THETA=5 deg 
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Figure 13: XG vs. K*GAMMA for THETA=0 deg 
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Figure 14: XG vs. K*GAMMA for THETA=-5 deg 
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Figure 15: XG vs. K*GAMMA for THETA=-10 deg 
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V. CONCLUSIONS AND RECOMMENDATIONS 


This thesis presented a continuing study of the formation of limit cycles due to the 
coupling of the sway/yaw/roll equations of motion. We have shown that loss of stability 
occurs in the form of stable limit cycles and that the addition of a nonzero pitch angle does 
not significantly affect the formation of these limit cycles. Through a linear analysis of the 
sway/yaw/roll equations of motion we demonstrated that the addition of a nonzero pitch 
angle affects the domain of stability of straight line motion, especially at higher speeds. 
This was validated by the nonlinear analysis as well. As a recommendation for further 
study in this area we suggest that the analysis be expandéd to include coupling into the 


heave/pitch directions of motion as well as the effects automatic path control. 
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APPENDIX 


The following is a list of the MATLAB and FORTRAN programs used in this thesis. 
Complete printouts of the programs accompany this list. 
e THESIS1.M 
A MATLAB program for performing the linear analysis section of this thesis. 
e HOPF_INEW.FOR 


A FORTRAN program for performing the nonlinear analysis section of this thesis. 


S? 


% THESIS 1.M 
%o 
% LOSS OF STABILITY 


J] 2K 2K 2K 2K 2K 2K 2k OK 2k 2K 2K EK OK 2k KOK OK BK OR OK OB OK 2 OK OK 2 KB KK OK Ok OK OK OK OK OK 


clear 


aaa 

W=12000; 

IXX=1760; TY Y=9450; 
IZZ=10700; [IXZ=0; LX Y=0; 
PY Z=0: L=17.425: RHO=1,94. 
G=32.2; U=8.0; M=W/G; B=W; 
THETA=0; 
THETA=THETA*pi/180; 
OMEGA=U*tan(THET A); 


NDI=034Re1© Ir 2. 
J. DEFINE HYDRODYNAMIC COEFFICIENTS 


YPDOT=1.270e-04*ND 1*L/2; 
YVDOT=-5.550e-02*ND1*L; 
YRDOT=1.240e-03*ND1*L’2; 

~ YP=3.055e-03*ND1*L; 
YPOMEGA=0; 
YP=YP*U+YPOMEGA*OMEGA+M*OMEGA; 
YV=-9.310e-02*ND1; 
YVOMEGA=0; 
YV=YV+YVOMEGA*OMEGA; 
-YR=-5.940e-02*ND1*L; 
YROMEGA=0; 
YR=YR+YROMEGA*OMEGA; 


NPDOT=-3.370e-05*ND1*L‘3; 
NVDOT=1.240e-03*ND1*L‘2; 
NRDOT=-3.400e-03*ND1*L43; 
NP=-8.405e-04*ND1*L‘2; 
NPOMEGA=0; 
NP=NP+NPOMEGA* OMEGA; 
NV=-1.484e-02*ND1*L; 
NVOMEGA=0; 
NV=NV+NVOMEGA*OMEGA; 
NR=-1.640e-02*ND1*L‘2; 
NROMEGA=0; 
NR=NR+NROMEGA*OMEGA; 


—_ == 
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KPDOT=-1.01le-03*ND1*L‘43; 
KVDOT=1.27e-04*ND1*L‘2; 
KRDOT=-3.37e-05*ND1*L43; 
=-]1.10e-02*ND1*L‘2; 
KPOMEGA=0; 
KP=KP+KPOMEGA*OMEGA; 
KV=3.055e-03*ND1*L; 
KVOMEGA=0; 
KV=KV+K VOMEGA*OMEGA; 
KR=-8.41e-04*ND1*L‘2; 
KROMEGA=0; 
KR=KR+KROMEGA*OMEGA; 


flag=0; 
for XG=0:0.01:2, 


flag=flag+1; 


xg(flag)=XG; 

a=IX X-KPDOT; b=KP*U; e=K V*U; 
=KRDOT 1=Y P+: 7=MeY VDOT; kK=YV+U; 
1]=XG*M-YRDOT; m=U*(YR-M); o=NPDOT; 
p=NP*U; gq=-XG*W; r=XG*M-NVDOT; 
w=U*(NR-XG*M); x=NV*U; u=IZZ-NRDOT; 


al=-u*¥M*%2; 
a2=-u*M* YPDOT-u*M*K VDOT+M*o0*]+f*r*M; 
a3=a*] *y-a*]*r-u*K VDOT* YPDOT+K VDOT *o0*l+f*r* YPDOT-f*0*}j; 


b1l=(w*(M42))-+(r*U*(M42)); 
b2=-M*e*u-M*u*i+w*M* YPDOT+w*K VDOT*#M- o*m*M4p*l*M- 
f*x*M+r*M*U* YPDOT-0*j*M*U+r*KR*U*M: 

b3=-e*u* YPDOT+e*o*l-u*i*K VDOT+w*K VDOT*YPDOT- 
KVDOT*o*m+K VDOT *p*l-a*j*w-a*k*uta*l*x-+a*r*m-b*j*u-+b*]*r+f*r*i- 
f*x* YPDOT+0*k*f-f*p*j+r*KR*U* YPDOT; 


b2=b2+M*OMEGA*j*u-M*OMEGA*I*r; 


cl=-x*(M‘2)*U; 

c2=r*1* M*U-x*M*U* YPDOT-x* KR*U*M-+0*k*M*U-p*]*#M*U04)*u* W-1*r*W- 
q*]*M+w*i*M-m*p*M+e*w*M; 
c3=a*k*w-a*m*x+b*j*w+b*k*u-b*]*x-b*r*m-41r*1* KR*U- 

X*KR*U* YPDOT+0*k* KR*U-p*j*KR*U-q*]* K VDOT+w*1* KVDOT-m*p*K VDOT- 
e*u*it+e*w* YPDOT-e*o*m+e*p*l+f*q*j-f*x*i+f*p*k; 


—OXW" 


4] 


c2=c2-M*OMEGA*j*W-M*OMEGA*k*u-M*OMEGA*]*x+M*OMEGA*r*m; 


d2=q*j*M*U-x*1*M*U+p*k*M*U+q*m*M-j]*w* W-k*u*W+41*#x*W4r*m*W; 
d3=q*j*KR*U-x*1* KR*U0+p*k*KR*U-f*q*k+q*m*K VDOT-e*q*l+e*w*i-e*m*p- 
b*k*w+b*¥m*x; 


d2=d2+M*OMEGA*k*w-M*OMEGA*U*(KR-M)*x; 


e2=k*w*W-m*x* W-q*k*M*U; 
e3=e*q*m-q*tk* KR*U; 


£5=(-e2*(b142))+b1*c1 *d2: 
£4=((b2*cl+b1 *c2)*d2)+b1 *c1 *d3-al *(d2“2)-e3*(b1%2)-2*e2*b1 *b2; 
£3=-a2*(d2%2)-2*d3*d2*al -2*e3*b 1 *b2- 

e2*((b242)+2*b 1 *b3)+d2*(b3*c1+b2*c2+b1 *c3)4+d3*(b2*cl+b1 *c2); 
f2=-e3*((b2%2)+2*b 1 *b3)-2*%e2*b2*b3+d2*(b3*c2+b2*c3)+d3*(b3*cl +b2*c2+b 1 *c3)- 
a3*(d2%2)-2*d3*d2*a2-al *(d342); 
f1=b3*c3*d2+d3*(b3*c2-+b2*c3)-2*e3*b2*b3-e2*(b3%2)-2*d3*d2*a3-a2*(d3%2); 
f0=b3*c3*d3-a3*(d3%2)-e3*(b3%2); 


coef=[f5 f4-f3 f2 f1 f0]; 
ZG=roots(coef); 


tot(flag)=ZG(5, 1); 
end 
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COE Ore 


Ce 


PROGRAM HOPF_1NEW.FOR 
EVALUATION OF HOPF BIFURCATION FORMULAS USING THE SUBOFF 
SUBMARINE MODEL 


IMPLICIT DOUBLE PRECISION (A-H,O-Z) 

REAL*8 LY Y,M,YPDOT,YVDOT,ND1,YRDOT 
REAL*8 YP, YV,YR,NPDOT,NVDOT,NRDOT,NP,NV,NR 
REAL*8 GAMA,U,KVDOT,KRDOT,KPDOT,D 

REAL*8 E0,E1,E2,E3,E4,XG,ZG,KR,KP,KV,XG1,ZG1 


REAL*8 M11,M12,M13,M14,M21,M22,M23 

REAL*8 M24,M31,M32,M33,M34,M41,M42,M43,M44 
REAL*8 N11,N12,N13,N14,N21,N22,N23,N24 
REAL*8 N31,N32,N33,N34,N41,N42,N43,N44 
REAR Ss LVI iesrieribeel 2 210235,024 031 

Re Awe ao leazemomriba4 Ios 161e1 7 125 1e26.077 235 
REAL*8 L36,L37,L1A,L2A,L3A,L4A,L5A,L6A,L7A,L8A,L9A 
REAL*8 L10A,011A,L12A,01,L2,L3,L4,L5,L6,L7 
REAL*8 L8,L9,R11,R12,R13,R14,R21,R22,R23,R24 
REA aS EAerI2Z-P13,P21,P22,P23 

DOUBLE PRECISION THETA 


DIMENSION F(4,4),T(4,4), TINV(4,4),FV 1(4),IV1(4), YYY(4,4) 

DIMENSION WR(4),W1(4), TSA VE(4,4), TLUD(4,4), IVLUD(4),SVLUD(4) 
DIMENSION ASAVE(4,4),A2(4,4),XL(18),HT(18),ZGG(202),FF(4,4) 
DIMENSION VEC0(18), VEC1(18), VEC2(18), VEC3(18), VEC4(18),XGG(202) 


INTEGER [J,K 


OPEN (20,FILE="HOPF25.RES') | 
OPEN (21,FILE='DAT25.DAT’ STATUS='OLD') 
OPEN (23,FILE="HOPF1.RES',STATUS='OLD') 


WEIGHT=12000.0 
IXX =1760.0 
TYY =9450.0 
IZZ =10700.0 
IXZ =0.0 

IXY =0.0 

YYZ =0.0 

L =17.425 
RHO =1.94 

CD =0.5 
CDOS eDaO.?)? 
G =32.2 


—_—i—=i< 
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XB =0.0 
ZB =0.0 
YG =0.0 
YB =0.0 
YDELTAR=0.0 
DELTAR=0.0 
NDELTAR=0.0 
NPROP=0.0 
M = WEIGHT/G 
B = WEIGHT 
W=WEIGHT 
WRITE (*,*) (ENPERSG! 
READ (*,*) U 
WRITE (*,*) ' ENTER THETA (DEGREES)' 
READ G} ) THETA 
U =5.0 
ND1 =0.5*RHO*L**2 
THETA=5 
THETA=THETA*PI/180 
OMEGA=U*DTAN(THETA) 


DEFINE HYDRODYNAMIC COEFFICIENTS 
YPDOT=1.270E-04*ND1*L**2 
Y¥ VDOT=-5.550E-02*ND1*L 
YRDOT=1.240E-03*ND1*L**2 
YP=3.055E-03*ND1*L 
YPOMEGA=0.D0 
YP=YP+YPOMEGA*OMEGA+M*OMEGA 
Y V=-9.310E-02*ND1 
YVOMEGA=0.D0 
YV=YV+Y VOMEGA*OMEGA 
YR=-5.940E-02*ND1*L 
YROMEGA=0.D0 
YR=YR+YROMEGA*OMEGA 


NPDOT=-3.370E-05*ND 1*L**3 
NVDOT=1.240E-03*ND1*L**2 
NRDOT=-3.400E-03*ND1*L**3 
NP=-8.405E-04*ND1*L**2 
NPOMEGA=0.D0 
NP=NP+NPOMEGA*OMEGA 
NV=-1.484E-02*ND1*L 
NVOMEGAS=0.D0 

~ NV=NV+NVOMEGA*OMEGA 
NR=-1.640E-02*ND1*L**2 


— ~~ 
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oO ene 


NROMEGA=0.D0 
NR=NR+NROMEGA*OMEGA 


KEDOt—-F0te-O3-NDIZL**3 
KVDOT=1.27E-04* NDI*L**2 
KRDOT=-3.37E-05*ND1*L**3 
Ke—lenOb-02* ND iss? 
KPOMEGA=0.D0 
KP=KP+KPOMEGA*OMEGA 
KV=3.055E-03*ND1*L 
KVOMEGA=0.D0 

KV=KV+K VOMEGA*OMEGA 
KR=-8.41E-04*ND 1*L**2 
KROMEGA=0.D0 
KR=KR+KROMEGA*OMEGA 


DEFINE THE LENGTH X AND HEIGHT H TERMS FOR THE INTEGRATION 
ALL IN FEET. 

XL( 1)=-105.9/12.0 
XL( 2)=-104.3/12.0 
XL( 3)=-99.3/12.0 
XL( 4)=-94.3/12.0 
XL( 5)=-87.3/12.0 
XL( 6)=-76.8/12.0 
XL( 7)=-66.3/12.0 
XL( 8)=-55.8/12.0 
XL( 9)=72.7/12.0 
XL(10)=79.2/12.0 
XL(11)=83.2/12.0 
XL(12)=87.2/12.0 
XL(13)=91.2/12.0 
XL(14)=95.2/12.0 
XL(15)=99.2/12.0 
XL(16)=101.2/12.0 
XL(17)=102.1/12.0 
XL(18)=103.2/12.0 


HT( 1)= 0.000 

HT( 2)= 2.28/12.0 
HT( 3)= 8.24/12.0 
HT( 4)= 13.96/12.0 
HT( 5)= 19.76/12.0 
HT( 6)= 25.1/12.0 
HT( 7)= 29.36/12.0 
HT( 8)= 31.85/12.0 


——~ 
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HT( 9)= 31.85/12.0 
HT(10)= 30.00/12.0 
HT(11)= 27.84/12.0 
HT(12)= 25.12/12.0 
HT(13)= 21.44/12.0 
HT(14)= 17.12/12.0 
HT(15)= 12.0/12.0 
HT(16)= 9.12/12 
HT(17)= 6.72/12 
HT(18)= 0.00 


DO 104 K = 1,18 
VEC0(K)=HT(K) 
VEC1(K)=XL(K)*HT(K) 
VEC2(K)=XL(K)*XL(K)*HT(K) 
VEC3(K)=XL(K)*XL(K)*XL(K)*HT(K) 
VEC4(K)=XL(K)*XL(K)*XL(K)*XL(K)*HT(K) 

104 CONTINUE 

CALL TRAP(18, VECO,XL,E0) 

CALL TRAP(18,VEC1,XL,E1) 

CALL TRAP(18, VEC2,XL,E2) 

CALL TRAP(18,VEC3,XL,E3) 

CALL TRAP(18,VEC4,XL,E4) 


GAMA=0.001 


— SS ee ee 
= ee 


C READ THE CRITICAL VALUES FOR XG AND ZG FROM FILE DATAI.DAT 


XGG(1)=0.0 

ZGG(1)=0.016358083 
DO 1 [=2,202 

READ (21,*)XG,ZG 


c WRITE(*,*)XG,ZG 
XGG(D=XG 
ZGG(1)=ZG 
CaSSSSSSS See Qa SSS SSS SSS SSS SSS SSS SR SS SSS SSS SSS S SSS SSSSSSSS SSS — = 
c DETERMINE [F] COEFFICIENTS 
C 


D=((M-Y VDOT)*(IZZ-NRDOT)*(IXX-KPDOT)) 

& -((M-Y VDOT)*(IXZ-KRDOT)*(-IXZ-NPDOT)) 

& -((M*XG-NVDOT)*(M*XG-YRDOT)*(IXX-KPDOT)) 
& +((M*XG-NVDOT)*(IXZ-KRDOT)*(-M*#ZG-YPDOT)) 
& +((-M*ZG-KVDOT)*(M*XG-YRDOT)*(-IXZ-NPDOT)) 


——~ 
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& -((-M*ZG-K VDOT)*(IZZ-NRDOT)*(-M*ZG-YPDOT)) 


A11=((IZZ-NRDOT)*(IXX-KPDOT))-((IXZ-KRDOT)*(-IXZ-NPDOT)) 
A12=((-M*XG+YRDOT)*(IXX-KPDOT))+((IXZ-KRDOT)*(-M#ZG-YPDOT)) 
A13=((M*XG-YRDOT)*(-IXZ-NPDOT))-((IZZ-NRDOT)*(-M*ZG-YPDOT)) 
A21=((-M*XG+NVDOT)*(IXX-KPDOT))+((-M*ZG-K VDOT)*(-IXZ-NPDOT)) 
A22=((M-YVDOT)*(IXX-KPDOT))-((-M*ZG-K VDOT)*(-M*ZG-YPDOT)) 
A23=((-M+Y VDOT)*(-IXZ-NPDOT))+((M*XG-NVDOT)*(-M*ZG-YPDOT)) 
A31=((M*XG-NVDOT)*(IXZ-KRDOT))-((-M*ZG-KVDOT)*(IZZ-NRDOT)) 
A32=((-M+YVDOT)*(IXZ-KRDOT))+((-M*ZG-K VDOT)*(M*#XG-YRDOT)) 
A33=((M-YVDOT)*(IZZ-NRDOT))-((M*XG-NVDOT)*(M*#XG-YRDOT)) 


EVALUATE TRANSFORMATION MATRIX OF EIGENVECTORS 


F(1,1)=(A11* YV*U+A i12*NV*U+A 13*KV*U)/D 
F(1,2)=(A11*(YR*U-M*U)+A12*(-M*XG*U+NR*U)+A 13*(M*ZG*U+KR*U))/D 
F(1,3)=(A11* YP*U+A 12*NP*U+A 13*(KP*U-M*ZG*OMEGA))/D 
F(1,4)=(A11*(W-B)*DCOS(THETA)+A12*XG*(W-XB)*B*DCOS(THETA)+ 

& A13*(-ZG*W+ZB)*B*DCOS(THETA))/D 

F(2,1)=(A21* YV*U+A22*NV*U+A23*KV*U)/D 
F(2,2)=(A21*(YR*U-M*U)+A22*(-M*XG*U+NR*U)+A23*(M*ZG*U+KR*U))/D 
F(2,3)=(A21* YP*U+A22*NP*U+A23*(KP*U-M*ZG*OMEGA))/D 
F(2,4)=(A21*(W-B)*DCOS(THETA)+A22*(XG*W-XB*B)*DCOS(THETA)+ 

& A23*(-ZG*W+ZB*B)*DCOS(THETA))/D 
F(3,1)=(A31* YV*U+A 32*NV*U+A33*KV*U)/D 
F(3,2)=(A31*(YR*U-M*U)+A32*(-M*XG*U+NR*U)+A33*(M*ZG*U+KR*U))/D 
F(3,3)=(A31* YP*U+A32*NP*U+A33*(KP*U-M*ZG*OMEGA))/D 
F(3,4)=(A31*(W-B)*DCOS(THETA)+A32*(XG*W-XB*B)*DCOS(THETA)+ 

& A33*(-ZG*W+ZB*B)*DCOS(THETA))/D 

F(4,1)=0.0 

F(4,2)=0.0 

F(4,3)=1.0 

F(4,4)=0.0 


DO 11 K=1,4 
DO 12 J=1,4 
ASAVE(K,J)=F(K,J) 


Wee CONTINUE 


11 CONTINUE 


CALL RG(4,4,F,WR,WLI1,YYY,IV1,FV1,IERR) 
WRITE(23, 1007) WR(1), WR(2), WR(3), WR(4) 
CALL DSOMEG(IEV,WR,WILOMEGA,CHECK) 
WRITE (*,*) CHECK 
WRITE (*,*) WR(2) 


—_—asa 
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C WRITE (*,*) WI(4) 
OMEGA0=OMEGA 
DO 5 J=1,4 
T(J,1)= YYY(,IEV) 
T(J,2)=-YYY (J IEV+1) 

5 CONTINUE 
IF (IEV.EQ.1.0) GO TO 13 
IF (IEV.EQ.2.0) GO TO 18 
IF (IEV.EQ.3.0) GO TO 14 

C STOP 3004 

14 DO6J=1,4 
T(J,3)=YYY(,1) 
T(,4)=YYY(J,2) 

6 CONTINUE 
GOm0 17 

18 DO 19 J=1,4 
T(I,3)=YYY(,1) 
T(LM=YYY(,4) 

19 CONTINUE 
Go unouy 

13 DO 16 J=1,4 
T,3)=YYY(J,3) 
T(I,4)=YYY(I,4) 

16 CONTINUE 

17 CONTINUE 


NORMALIZATION OF THE CRITICAL EIGENVECTOR 


) GO 


CALL NORMAL(T) 


INVERT TRANSFORMATION MATRIX 


Gee 


DO 2 K=1,4 
DO 3 J=1,4 
TINV(K,]J)=0.0 
TSAVE(K,J)=T(K,J) 
3. CONTINUE 
2 CONTINUE 
CALL DLUD(4,4,TSAVE,4,TLUD,IVLUD) 
C DpO4J=1,4 
Ee IF (IVLUD(J).EQ.0) STOP 3003 
C 4 CONTINUE 
CALL DILU(4,4, TLUD,IVLUD,SVLUD) 
DO 8 K=1,4 
DO 9 J=1,4 
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TINV(K,J)=TLUD(K,J) 


9 CONTINUE 
8 CONTINUE 
C 
C CHECK Inv(T)*A*T 
C 
CALL MULT(TINV,ASAVE,T,A2) 
@ 
P1=A2(3,3) 
P2=A2(4,4) 
PEIG1=P1 
PEIG2=P2 
e 
C DEFINITION OF Nij 
GE 


N11=TINV(1,1) 
N1Z=TINV(1,2) 
N13=TINV(1,3) 
N14=TINV(1,4) 
N21=TINV(2,1) 
N22=TINV(2,2) 
N23=TINV(2,3) 
N24=TINV(2,4) 
N31=TINV(3,1) 
N32=TINV(3,2) 
N33=TINV(3,3) 
N34=TINV(3,4) 
N41=TINV(4,1) 
N42=TINV(4,2) 
N43=TINV(4,3) 
N44=TINV(4,4) 


DEFINITION OF My 


OFGiee 


M11=T(1,1) 
M12=T(1,2) 
M13=T(1,3) 
M14=T(1,4) 
M21=T(2,1) 
M22=T(2,2) 
M23=T(2,3) 
M24=T(2,4) 
M31=T(3,1) 
M32=T(3,2) 
M33=T(3,3) 


——* 


M34=T(3,4) 
M41=T(4,1) 
M42=T(4,2) 
M43=T(4,3) 
M44=T(4,4) 


DEFINITION OF Liy 


Or) se) 


L1=YG*((M31**2)+(M21**2)) 

L2=IXY*(M31**2)4+1TYZ*M31*M21+YG*M21*M11 

L3=IX Y*M31*M214+1YZ*(M21**2)-+M*YG*M11*M31+YG*W*(M41**2)- 
YB*B*(M41** 

&2) | 

LA=YG*(2.0*M31*M32+2.0*M21*M22) 


L5=2.0* IX Y*M31*M32+TYZ* (M3 1*M224+M32*M21)+YG*(M21*M12+M22*M11) 


L6=LX Y *(M31*M22+M32*M21)+2.0*TYZ*M21*M22+M* YG*(M11*M32+M12*M31 
)+2. 

&0* (CY G* W- YB*B)*M41*M42 

72 GaAs 2**2)-4(M22**2)) 

L8=IX Y*(M32**2)+TYZ*M32*M22+YG*M22*M12 

L9O=IX Y*M32*M22+7Y Z*(M22**2)+M* YG*M12*M32+(YG*W- 
YB*B)*(M42**2) 
Cc 
Cc 
L15=(A11/D)*L1+(A12/D)*L2-(A13/D)*L3 
L16=(A11/D)*L4+(A12/D)*L5-(A13/D)*L6 
L17=(A11/D)*L7+(A12/D)*L8-(A13/D)*L9 
L25=(A21/D)*L1+(A22/D)*L2-(A23/D)*L3 
L26=(A21/D)*L4+(A22/D)*L5-(A23/D)*L6 
L27=(A2 1/D)*L7+(A22/D)*L8-(A23/D)*L9 
L35=(A31/D)*L1+(A32/D)*L2-(A33/D)*L3 
L36=(A31/D)*L4+(A32/D)*L5-(A33/D)*L6 
L37=(A31/D)*L7+(A32/D)*L8-(A33/D)*L9 


C=CD/(6.0*GAMA) 


L1A=C*(E0*(M11**3)+3.0*E 1 *(M11**2)*M214+3.0*E2*M 1 1*(M21**2)+E3*(M21 
&**3))-((1.0/6.0)*(W-B)*DCOS(THETA))*(M41*#3) 


L2A=C*(B1*(M11**3)+3.0*E2*(M11**2)*M21+3.0*E3*M1 1*(M21**2)+E4*(M21 
&**3))-((1.0/6.0)*(XG* W-XB*B)*DCOS(THETA))*(M41**3) 
L3A=((1.0/6.0)*(ZG* W-ZB*B)*DCOS(THETA))*(M41**3) 
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TAAK GeO 0 (M11 *2) = M2307 1 (Mil **2)*#M22--24M) 1*M12*M21)+3.0 
&*E2*(M12*(M21**2)+2.0*M21*M22*M1 1)+3.0*E3*(M21**2)*M22)-((1.0/6.0) 
&*(W-B)*DCOS(THETA))*3.0*(M41**2)*M42 


LSA=C*(3.0*E1*(M 11**2)*M12+3.0*E2*((M11**2)*M22+2*M11*M12*M21)+3.0 
& *E3*(M12*(M21**2)+2*M21*M22*M1 1)+3.0*E4*(M21**2)*M22)-((1.0/6.0)*(X 
&G*W-XB*B)*DCOS(THETA))*3.0*(M41**2)*M42 
L6A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*3.0*(M41**2)*M42 


L7A=C*(3*EO*M 1 1*(M12**2)+3.0*E1*((M12**2)*M21+2*M11*M12*M22)+3.0*E 
&2*((M22**2)*M11+2.0*M21*M22*M12)+3.0*E3*M21*(M22**2))-((1.0/6.0)*(W 
&-B)*DCOS(THETA))*3.0*M41*(M42**2) 


L8A=C*(3.0*E1*M11*(Mi2**2)+3.0*E2*((M12**2)*M21+2.0*M11*M12*M22)+3 
&.0*E3*((M22**2)*M114+2*M21*M22*M12)+3.0*E4*M21*(M22**2)) 
L9A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*3.0*M41*(M42**2) 
L10A=-C*(E0*(M12**3)+3.0*E1*(M11**2)*M214+3.0*E2*M 12*(M22**2)+E3*(M 
&22**3))-((1.0/6.0)*(W-B)*DCOS(THETA))*(M42**3) 
L11A=-C*(E1*(M12**3)+3.0*E2*(M11**2)*M21+3.0*E3*M 12*(M22**2)+E4*(M 
&22**3))-((1.0/6.0)*(XG*W-XB*B)*DCOS(THETA))*(M42**3) 
L12A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*(M42**3) 


L11=(-A11/D)*L1A+(A12/D)*L2A+(A13/D)*L3A 
L12=(-A11/D)*L4A+(A12/D)*L5A+(A13/D)*L6A 
L13=(-A11/D)*L7A+(A12/D)*L8A+(A 13/D)*L9A 
L14=(-A11/D)*L10A+(A12/D)*L11A+(A13/D)*L12A 


L21=(-A21/D)*L1A+(A22/D)*L2A+(A23/D)*L3A 
L22=(-A21/D)*L4A+(A22/D)*L5A+(A23/D)*L6A 
L23=(-A21/D)*L7A+(A22/D)*L8A4+(A23/D)*L9A 
L24=(-A21/D)*L10A+(A22/D)*L11A+(A23/D)*L12A 


L31=(-A31/D)*L1A+(A32/D)*L2A+(A33/D)*L3A 
L32=(-A31/D)*L4A+(A32/D)*L5A+(A33/D)*L6A 
L33=(-A31/D)*L7A+(A32/D)*L8A+(A33/D)*L9OA 
L34=(-A31/D)*L10A+(A32/D)*L11A+(A33/D)*L12A 


R11=(N11*L11)+(N12*L21)+(N13*L31) 
R12=(N11*L12)+(N12*L22)+(N13*L32) 
R13=(N11*L13)+(N12*L23)+(N13*L33) 
R14=(N11*L14)+(N12*L24)+(N13*L34) 
R21=(N21*L11)+(N22*L21)+(N23*L31) 
R22=(N21*L12)+(N22*L22)+(N23*L32) 
R23=(N21*L13)+(N22*L23)+(N23*L33) 
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@ 0) ©) 


Oe) 


ee) a: 


R24=(N21*L14)+(N22*L24)+(N23*L34) 


P11=(N11*L15)+(N12*L25)+(N13*L35) 
P12=(N11*L16)+(N12*L26)+(N13*L36) 
P13=(N11*L17)+(N12*L27)+(N13*L37) 
P21=(N21*L15)+(N22*L25)+(N23*L35) 
P22=(N21*L16)+(N22*L26)+(N23*L36) 
P23=(N21*L17)+(N22*L27)+(N23*L37) 


EVALUATE DALPHA AND DOMEGA 


ZGR =ZGG(D 
ZGL =ZGG(I-1) 
ZG1 =ZGR 
XG1 =XGG(I) 


CALL FMATRIX(ZG1,XG1,FF,U,THETA) 


CALL RG(4,4,FF,WR,W10,YYY,IV1,FV1,IERR) 
CALL DSTABL(DEOS,WR,WLFREQ) 
ALPHR=DEOS 

OMEGR=FREQ 


ZG1 =ZGL 
XG1 =XGG(I-1) 


CALL FMATRIX(ZG1,XG1,FF,U,THETA) 


CALL RG(4,4,FF,WR,WI10,YYY,IV1,FV1,IERR) 
CALL DSTABL(DEOS,WR,WLFREQ) 
ALPHL=DEOS 

OMEGL=FREQ 


DALPHA=(ALPHR-ALPHL)/(ZGR-ZGL) 
DOMEGA=(OMEGR-OMEGL)/(ZGR-ZGL) 


EVALUATION OF HOPF BIFURCATION COEFFICIENTS 


COEF 1=(1.0/8.0)*(3.0*R11+R13+R22+3.0*R24) 
COEF2=(1.0/8.0)*(3.0*R11+R23-R12-3.0*R14) 


AMU?2 =-COEF1/(8.0*DALPHA)  ??72772? 
BETA2=0.25*COEF1 Re 

TAU2 =-(COEF2-DOMEGA*COEF1/DALPHA)/(8.0* OMEGAQO) 
PER =2.0*3.1415927/OMEGAO 


— O—~ 


a2 


C PER =PER*U/L 
WRITE (20,2001)XG,ZG,COEF 1, DALPHA,OMEGAO,PEIG1,PEIG2 
C WRITE (20,2001)COEF1 
1 CONTINUE 
c 
STOP 


2001 FORMAT (7E14.5) 

4001 FORMAT (3F15.5) 

1007 FORMAT (4E14.5) 
END 


5S 
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